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ABSTRACT 

''J 

We  the  problem  of  boundary  conditions  for  the  vortidty  form  of  the  2-D 

Navier-Stokes  equations  and  the  Prandtl  boundary  layer  equations.  We  present  a  new  formu¬ 
lation  of  the  vortidty  boundary  conditions  and  relate  these  conditions  to  those  which  are 
currently  used  in  mwwwncal  algorithms.  We  present  a  finite  difference  scheme  incorporating 
the  new  boundary  conditions  for  the  Prandtl  boundary  layer  equations.  Numerical  results  are 
presented. 
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[  1.  Vortidty  Boundary  CondMn—  for  tbo  2*D  Navier-Stokes  Equations 

i 

|  In  this  paper  we  present  some  preliminary  results  concerning  investigations  into  the  issue 

of  boundary  conditions  for  vortidty  and  the  boundary  generation  of  vortidty  in  the  context  of 
the  two  dimensional  incompressible  Navier-Stokes  equations  and  the  Prandtl  boundary  layer 
equations.  We  first  discuss  the  origin  of  the  problem  with  vortidty  boundary  values  for  the 
Navier-Stokes  equations.  We  then  discuss  two  techniques,  one  due  to  Chorin  [3],  and  the 
other  due  to  QuartapeUe  and  Valz-Gris  [8],  for  overcoming  this  problem.  These  two  tech* 

i 

I  niques  seem  very  different,  and  yet' when  incorporated  into  numerical  schemes,  both  appear  to 

give  reasonable  solutions  to  the  Navier  Stokes  equations.  One  of  the  primary  motives  for  this 
investigation  was  to  understand  why  this  should  be  so.  In  order  to  obtain  this  understanding, 
we  found  it  necessary  to  derive  our  own  solution  to  the  problem  of  boundary  conditions.  We 
present  this,  and  use  the  derivation  to  help  explain  the  relationship  between  the  two  earlier 
approaches.  Our  derivation  also  suggests  new  methods  for  implementing  the  vortidty  boun¬ 
dary  conditions  in  numerical  schemes.  This  aspect  will  be  discussed  in  a  future  paper,  how¬ 
ever,  to  illustrate  the  general  idea  of  our  approach  we  consider  an  application  of  our  technique 
to  the  Prandtl  boundary  layer  equations.  We  present  a  numerical  scheme  for  these  equations 
and  provide  some  computational  results  obtained  with  this  scheme.  We  also  compare  a  finite 
difference  scheme  for  these  equations  with  a  scheme  introduced  by  Chorin  in  [4]. 

The  equations  we  are  concerned  with  are  the  two  dimensional  incompressible  Navier- 
Stokes  equations, 

4^  +  a  ■  va  =  -vp  +  vAu  (l.i) 

of 

v  ■  a  •  o  (1.2) 

tf(a)  =  8(a)  for  a  c  <30  .  (13) 

Here  il  is  the  velotity,  P  the  pressure  and  v  the  viscosity.  We  assume  the  fluid  is  of  constant 
density  equal  to  one.  O  is  the  region  in  R2  with  boundary  <30.  8(a)  is  the  veloaty  on  the 
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boundary,  often  taken  to  be  identically  zero. 

In  the  vartidty  formulation  of  (1.1)-(1.3)  the  velocity  field  a  is  taken  to  be  the  sum  of  a 
velocity  field  due  to  an  irrotational  flow  and  a  velocity  field  due  to  a  rotational  flow.  Let  <J>  be 
the  potential  for  the  irrotational  flow  and  Y  die  stream  function  for  the  rotational  flow.  We 
assume  4>  satisfies 

■  0  *  8  •  A  on  9(1  . 

dii 

Here  A  denotes  the  normal  to  the  boundary 

Let  cd  be  the  vortidty.  By  taking  die  curl  of  equation  (1.1)  one  obtains  the  following 
equation  for  the  transport  of  vortidty, 

+  U  •  Vw  *  vAw  .  (1.4) 

tT - 

not 

u  -  (<►„  ♦,)  +  (%,  -Yt)  (1.5) 

and  Y  is  determined  from  the  equation 

AY  *  -w 

Y  *  0  on  dfl  .  (1.6) 

The  boundary  conditions  given  above  for  $  and  Y  guarantee  that  the  normal  veloaty 
boundary  condition  is  To  satisfy  the  tangential  veloaty  condition  we  must  also  have 

~  *  S  •  f  -  (d»T,4»v)  •  f  =  b(x)  for  x  e  aft  (1.7) 

where  f  is  the  unit  tangent  to  dfl. 

Often  the  vortidty  form  of  the  equations  is  used  in  numerical  computations.  One  reason 
for  its  use  is  that  the  number  of  unknowns  is  reduced  from  three  to  one.  Also,  in  many  flows 


the  vorticity  is  itself  of  great  interest  and  it  is  desirable  to  compute  it’s  evolution  directly 
rather  than  compute  in  the  primitive  variables  and  differentiate  the  result.  The  use  of  the  vor¬ 
tidty  formulation  is  not  without  it’s  difficulties.  One  major  difficulty  is  the  proper  vortidty 
boundary  conditions. 

In  the  translation  of  the  equations  from  primitive  variables  (1.1)-(1.3)  to  the  vorticity 
form  (1.4)-(1.7),  boundary  conditions  for  the  vorticity  are  not  obtained.  Moreover,  it  appears 
that  too  many  boundary  conditions  are  given  on  the  stream  function.  (Both  those  in  (1.6)  and 
(1.7)  must  be  satisfied.)  However,  a  bit  of  thought  reveals  that  the  freedom  in  choosing  the 
vorticity  boundary  conditions  occurs  precisely  because  the  stream  function  is  over  determined. 
For,  unless  there  is  some  mechanism  for  manipulating  u  in  the  interior  of  the  domain,  the 
problem  which  Hrfiw  qr  will  not,  in  general,  be  well  posed.  It  is  therefore  not  surprising 
that  a  common  thread  which  runs  through  all  numerical  implementations  is  the  manipulation 
of  the  vortidty,  usually  near  the  boundary,  in  a  way  which  insures  that  both  (1.6)  and  (1.7) 
are  satisfied.  However,  apart  from  this  similarity,  the  numerical  techniques  for  overcoming 
this  difficulty  can  appear  quite  different.  We  shall  consider  two  approaches  in  use.  There  are 
others,  most  notably  those  which  use  the  relation  between  Y  and  u>  in  (1.6)  and  (1.7)  at 
points  on  the  boundary.  For  a  review  of  these,  and  of  the  results  which  can  be  obtained  with 
them,  see  [6]  and  [7]. 

One  approach  to  solving  the  problem  of  boundary  conditions  is  that  which  is  implidtly 
used  by  Qwrin  in  [3].  In  that  paper,  he  utilizes  the  method  of  fractional  steps  to  solve  equa¬ 
tions  (1.4)-(1.7).  The  method  is  to  advance  for  one  timestep  the  solution  to  the  Euler  equa¬ 
tions: 

^  +  0  •  Vca  =  0  (1.8) 

(with  U  defined  by  equation  (1.5)  and  ¥  determined  only  by  (1.6))  followed  by  one  timestep 
of  the  solution  of  the  heat  equation 
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7T  “  •  (1-9) 

Die  vortex  blob  method  is  used  to  advance  the  solution  for  the  Euler  step,  and  a  random  walk 
technique  is  used  to  simulate  the  diffusion  step.  In  the  solution  of  the  Euler  equations,  (1.8), 
boundary  conditions  on  the  vortidty  are  not  needed  since  U  •  if  =  0.  In  the  solution  of  the 
heat  equation  (1.9)  boundary  conditions  are  needed,  and  one  of  the  novel  ideas  contained  in 
Chorin’s  paper  is  the  method  by  which  the  boundary  values  of  vortkity  are  obtained. 

After  one  step  of  the  Euler  equations,  the  interior  vortidty  distribution  will  not  neces¬ 
sarily  induce  a  velocity  field  which  satisfies  the  tangential  velocity  boundary  condition.  If  u  is 
the  distribution  of  vortidty  after  one  time  step  of  Euler,  and  (a,  v)  is  the  velodty  constructed 
from  this  distribution,  then  (£,  v)  •  f  #  b(s)  for  s  c  dfl.  Chorin’s  idea  is  to  introduce  vorti¬ 
dty  on  the  boundary  in  an  amount  which  cancels  this  error.  Specifically,  if  the  boundary  is 
divided  up  into  segments  of  length  A  then  die  vortidty  introduced  per  segment  is 

<•>(*/)  -  (a,  v)|1(  •  f  h  -  b(s,) 

where  r,  is  the  center  of  the  ith  segment  on  the  boundary.  This  vortidty  is  then  allowed  to  dif¬ 
fuse  into  the  region  and  participate  in  the  evolution  of  the  flow. 

A  common  concept  in  fluid  mechanics  is  that  vortidty  is  "produced”  at  object  boun¬ 
daries.  As  discussed  by  Batchelor  in  [2],  this  production  occurs  so  that  the  external  or  outer 
fluid  velodty  will  match  the  velodty  of  the  object.  An  attractive  feature  of  Chorin’s  method 
is  that  it  mimics  this  process  of  boundary  generation  of  vortidty.  Computations  performed 
with  this  method  (and  an  improved  version  incorporating  boundary  layer  modelling  [4])  sug¬ 
gest  that  this  algorithm  generates  a  good  approximation  to  the  solution  of  the  equations 
(1.4)-(1.7).  From  this  fact  we  conclude  that  boundary  vortidty  generation,  as  understood 
from  the  nature  of  the  physical  problem,  is  intimately  related  to  the  missing  boundary  condi¬ 
tions  for  the  vortidty  transport  equations. 
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The  second  approach  is  that  due  to  QuartapeUe  and  Valz-Gris  [8]  and  is  more  mathemat¬ 
ically  rather  than  physically  motivated.  Essentially,  they  pose  the  following  question  "What 
are  the  conditions  on  the  vortidty  so  that  the  stream  function  constructed  using  (1.6)  also 
satisfies  (1.7)?”  The  answer  is  contained  in  the  following  theorem  presented  in  their  paper, 

Theorem.  A'P  -  -u>  in  ft,  ¥|,  =  a(s)  and - -  b(s),  if  and  only  if 

d* 

/<*>  r\da  =  /tih(r)  -  a(s)Qds  (1.10) 

n  »ci  dn 

for  any  function  q  harmonic  in  Cl. 

fit  particular,  if  a(s)  =  0  and  a  no  slip  condition  is  specified,  b(s)  =  0,  then  the  vorticity 
must  be  orthogonal  to  all  harmonic  functions  defined  in  the  domain. 

Using  this  result,  QuartapeUe  and  Valz-Gris  overcome  the  difficulty  of  finding  explicit 
boundary  conditions  on  the  vorticity  by  adjoining  the  constraint  (1.10)  to  the  equations  for 
vortidty  transport  (1.4).  Their  implementation  consists  of  advancing  the  solution  of  (1.4)  one 
time  step  without  regard  to  the  constraint  and  then  projecting  the  result  onto  the  component 
which  satisfies  (1.10).  It  is  for  this  reason  that  we  refer  to  their  approach  as  the  projection 
method.  Details  of  their  numerical  implementation  and  some  nice  computational  results  con¬ 
cerning  the  driven  cavity  problem  are  presented  in  [9]. 

Aside  from  the  fact  that  both  approaches  are  used  to  solve  the  same  set  of  equations, 
the  techniques  do  not  appear  to  be  related.  Chorin’s  scheme  effects  the  vortidty  on  the  boun¬ 
dary,  while  in  the  method  of  QuartapeUe  and  Valz-Gris,  the  interior  vortidty  is  effected  as 
weU.  In  order  to  relate  the  two,  we  found  it  necessary  to  obtain  a  precise  set  of  boundary 
conditions  for  the  vortidty.  Our  derivation  is  motivated  by  the  observation  of  QuartapeUe 
and  Valz-Gris  that  one  should  consider  the  evolution  of  the  vortidty  as  a  constrained  evolu¬ 
tion.  Instead  of  explidtly  finding  the  constraint  on  the  vortidty  and  adjoining  it  to  the  equa¬ 
tions  as  QuartapeUe  and  Valz-Gris  do,  we  precede  by  finding  conditions  which  insure  that  as 
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the  vortidty  evolves  the  constraint  will  automatically  be  satisfied. 

We  begin  by  expressing  the  "extra"  boundary  condition  (1.7)  as  a  constraint  on  the  vorti¬ 
dty.  Far  this  purpose  we  use  the  Green’s  function,  G(xrr),  for  the  domain  Cl.  This  function 
is  the  solution  of 

AG(xrf)  *  8(j)  x  c  ft  GCxpt)  =  0  x  c  dCl 
where  8(j)  is  a  Dirac  delta  function  located  at  s  c  A.  We  have 

TO  - 

n 

so  that  condition  (1.7)  can  be  written  as 

-f-J  G(x,i)«(jr,r)<&  -  -  TO  for  x  c  30  (1.11) 

an  n 

To  find  boundary  conditions  which  insure  that  solutions  of  (1.4)  induce  a  stream  func¬ 
tion  which  satisfies  both  (1.6)  and  (1.7),  we  find  boundary  conditions  which  guarantee  that  the 
derivative  with  respect  to  time  of  the  constraint  (1.11)  vanishes.  We  require 

a7(irjr  G(xM,,t)ds  +  TO)  -  0  for  x  « an 

If  we  use  the  fact  that  w  is  a  solution  of  (1.4)  and  the  Green’s  identity 

u  *  / G(xrj)Aw(j,t)dj  -  / -^■(*,a)<i>(a,r)da 
n  an dn 

we  find 

-  H  - 0 

-  G(xri)(-0  •  V«*»(x,r)  +  vA<u (s,t))ds  —  0 

dn  n 


-  V-J-/  G(xj)Ati>(j,t)ds  -  J-J  G(xj)U  •  Vta(s,t)ds 

an  n  °n  n 

-  ft  +  aM«  •  **4*  (1.12) 

Ibis  is  an  integral-differential  equation  which  determines  the  boundary  values.  Values 
satisfying  (1.12)  and  used  as  data  for  (1.4)  will  insure  that  the  time  derivative  of  (1.11)  will  be 
zero.  If  we  assume  that  the  initial  vortidty  satisfies  (1.11),  then  we  will  have  insured  that  this 
constraint,  and  hence  (1.7),  is  satisfied  for  all  time. 

If  one  uses  the  boundary  conditions  (1.12)  with  (1.4)-(1.6)  then  the  resultant  vortidty 
will  also  satisfy  (1.7),  and  hence  will  satisfy  the  constraint  as  formulated  by  Quartapelle  and 
Valz-Gris,  (1.10).  Conversely  if  one  solves  (1.4)  with  (1.10)  explicitly  satisfied  then  the  vorti¬ 
dty  will  satisfy  (1.12).  We  condude  that  Oiartapelle  and  Valz-Gris,  by  enforcing  (1.10) 
explidtly  at  every  timestep,  are  impliatly  implementing  the  boundary  conditions  (1.12). 

The  connection  between  the  boundary  conditions  (1.12)  with  the  technique  due  to  Qio- 
rin  is  less  dear.  If  one  considers  finite  difference  methods  for  (1.4)-(1.6)  using  (1.12),  then 
one  finds  that  they  have  the  structure  that  the  boundary  vortidty  is  obtained  by  solving  an 
integral  equation  with  a  forcing  function  which  is  proportional  to  the  error  in  the  tangential 
velodty  induced  by  one  step  of  Euler  flow.  Rather  than  get  into  the  details  of  the  numerical 
schemes,  this  result  can  by  illustrated  by  considering  an  explidt  scheme  which  is  discrete  in 
time  and  continuous  in  space: 

- — r — —  =  -U  •  Vu*  +  vA w*  (1.13) 

of 

Boundary  values  of  vortidty  are  needed  in  the  evaluation  of  the  Lapladan  on  the  right 
hand  side  of  (1.13).  These  values  can  be  approximated  using  (1.12).  The  right  hand  side  of 
(1.12)  is  evaluated  using  the  solution  at  time  kit  and  then  the  left  hand  side  is  solved  for  the 
boundary  values.  If  one  uses  the  expression  for  one  step  of  Euler  flow 


it 

to  approximate  the  convective  terms  in  the  right  hand  side  of  (1.12),  one  finds 


Hrf  ^  °<"x 

Thus  the  right  hand  side  of  (1.12)  has  the  interpretation  of  being  the  difference  in  the  tangen¬ 
tial  velocity  at  the  boundary  between  that  at  time  JfcBr  and  that  induced  by  one  step  of  Euler 
flow.  If  the  tangential  velocity  boundary  condition  is  satisfied  at  time  kit,  then  the  right  hand 
side  is  just  proportional  to  the  error  in  velocity  induced  by  one  step  of  Euler  flow. 

The  incorporation  of  boundary  conditions  (1.12)  has  the  flavor  of  that  utilized  by  Cho- 
rin,  however,  in  Charin’s  scheme  one  does  not  solve  an  integral  equation  to  determine  the 
boundary  vortkity.  We  imagine  that  Chorin's  method  of  constructing  the  strength  is  an 
approximation  to  the  inversion  of  the  integral  equation  on  the  left  hand  side  of  (1.12).  A 
detailed  analysis  has  yet  to  be  carried  out.  For  a  related  set  of  equations,  the  Prandtl  boun¬ 
dary  layer  equations,  it  is  very  dear  that  Chorin’s  vortidty  creation  is  equivalent  to  the  imple¬ 
mentation  of  vortidty  boundary  conditions.  For  a  further  discussion  of  this  latter  case,  see 
the  following  section. 

>0  condude,  one  should  view  the  process  of  solving  the  Navier-Stokes  equations  as  a 
problem  of  solving  a  set  of  equations  subject  to  a  constraint.  QuartapeUe  and  Valz-Gris’s 
technique  is  to  expliatly  adjoin  the  constraint  to  the  vortidty  evolution  equations.  Their 
numerical  implementation  takes  the  form  of  a  projection  method.  Alternatively,  one  can 
incorporate  the  boundary  conditions  derived  above.  These  conditions  will  insure  that  if  the 
constraint  is  satisfied  initially,  then  it  will  be  satisfied  for  all  time.  We  believe  Chorin's  scheme 
to  be  an  approximation  to  these  boundary  conditions.  Thus  the  schemes  have  the  same  goal, 
that  of  evolving  the  vortidty  in  such  a  way  that  the  stream  function  satisfy  boundary 
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conditions  (1.6)  and  (1.7),  but  they  differ  in  the  manner  in  which  this  goal  is  achieved.  This 
observation  explains  why  such  different  techniques  applied  to  the  same  set  of  equations  can 
both  yield  reasonable  solutions. 
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2.  Results  for  the  Prandtl  Boundary  Layer  Equations 

In  this  section  we  apply  the  ideas  of  the  first  section  to  the  Prandtl  boundary  layer  equa¬ 
tions.  This  set  of  equations,  which  we  shall  use  to  describe  laminar  flow  over  a  half-infinite 
flat  plate,  is  a  simpler  set  of  equations  than  the  2-D  Navier  Stokes  equations,  and  yet  when 
expressed  in  the  vortidty  formulation  possesses  the  same  problems  with  vortidty  boundary 
conditions.  We  shall  derive  the  continuous  boundary  conditions  for  the  vortidty,  and  then 
implement  a  numerical  method  using  these  boundary  conditions.  Results  of  computations  with 
this  mimwial  method  will  be  presented.  A  similar  discussion  concerning  the  implementation 
of  numerical  methods  for  the  2-0  Navier  Stokes  equations  will  be  presented  in  a  forth  coming 
paper  [1]. 

hi  the  vortidty  form,  the  Prandtl  boundary  equations  are 


*2L  +  a.Vtamv£2L 
at  ay2 

(2.1) 

i 

i 

3 

(2.2) 

S3 

II  *  U0  +  /ci)(jr j,t)ds 

y 

(2.3) 

v ,  -}**&*£* 

0  ox 

(2.4) 

with  boundary  conditions 


«  =  v  ■  0  at  y  3  0  x  >  0 


u  =  UD  at  y  =  +« 


(2.5) 

(2.6) 


Here  (u,v)  is  the  velodty  and  u  the  vortidty.  The  domain  is  the  quarter  plane  Osis« 
and  0  sjis», 


11 


As  is  the  case  with  the  Navier  Stokes  equations,  when  one  reconstructs  the  velocity  field 
from  the  vortidty,  it  is  not  automatic  that  the  boundary  condition  on  u  in  (2.5)  will  be  satis¬ 
fied.  Similarly,  boundary  conditions  necessary  to  dose  equation  (2.1)  are  not  given.  We  pre¬ 
cede,  as  in  the  first  section,  to  derive  boundary  conditions  by  expressing  the  condition  on  u  in 
(2.5)  as  a  constraint  on  the  vortidty  and  then  setting  time  derivative  of  this  constraint  equal  to 
zero.  Using  (2.3),  the  constraint  on  the  vortidty  is 


U0  +  J <a(xj,t)ds  =0  x  >  0  and  t  as  0 


Thus  we  require 


j-fm(xj,t)ds  =  0 
at  q 

-  f~(xj,t)dt  *  0 
o  <*• 

-  f—U  •  V«  +  v-— -- ds  —  0 
o  By2 

w 

-  HrUo  -  ■  (2.8) 

By  •  v-J 

(2.8)  is  the  desired  boundary  condition.  We  now  consider  a  numerical  method  for  solving 
these  equations  which  incorporates  these  boundary  conditions. 

Our  computational  domain  is  the  rectangular  region  described  by  the  points  (x,  v)  such 
that  0  s  is  and  0  s  y  s  y„.  The  mesh  we  use  is  rectangular  with  widths  dx  and  dy  in 
the  x  and  y  direction  respectively.  The  values  of  the  vortidty  and  velodty  are  computed  at 
the  grid  points  (idx,  jdy )  and  are  designated  by  c o,j  utJ  and  vtJ  with  0  s  i  s  m  and 
Osjs«.  To  approximate  (2.1)  we  use  a  one  step  expliat  method  (Euler’s  method)  to 
advance  the  solution  in  time,  and  approximate  the  advection  term  d  ■  Vw  using  a  second  order 
upwind  differencing  scheme  due  to  Colella  [5].  (We  believe  the  results  are  relatively 


■*  V.'Aw 
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independent  of  the  advection  scheme  used)  The  second  derivative  term  was  approximated  by 
central  differences.  Our  scheme  is  thus  expressed  by 

.  A,/**.  «*)  +  d;d r»tj 

where  A^(il*,  w*)  is  the  second  order  approximation  to  the  convective  term  in  (2.1)  and 


D'D-»f,  -  -  H  -1-  SL1 

>  >  ,J  dy2 


(2.9) 


To  construct  u  we  use  the  trapezoidal  rule  to  approximate  the  integral  in  (2.3).  We 
assume  no  vortidty  above  the  line  y  =  y„,  and  hence  u  *  U0  fat  pants  above  this  line.  For 
the  remaining  points  (j  as  n)  we  use 


UU 


-uo+  i 

p*J  2 


(2.10) 


v  is  approximated  using 


M  ,  »  _  V  l(Ui+U>  ~Ut-U>  +  "I*  l**l  "  ft 

iJ  7- i2'  2dx  '  7 


2 dx 


(2-11) 


i.e.  a  trapezoidal  rule  approximation  to  (2.4)  in  which  central  differences  are  used  to  approxi¬ 
mate  the  derivatives  of  u.  At  points  on  the  left  and  right  computational  boundaries  second 
order  one  sided  differences  are  used  in  (2.11)  instead  of  central  differences. 


It  remains  to  specify  the  boundary  conditions  used  for  the  vortidty.  At  the  left  edge  of 
the  computational  boundary,  the  inflow  side,  we  assumed  that  there  is  no  vortidty  immedi¬ 
ately  upstream.  At  the  right  edge  of  the  computational  boundary  we  assumed  outflow  boun¬ 
dary  conditions  -  i.e.  the  vortidty  is  unspecified  at  pants  just  outside  of  the  computational 
domain.  At  the  top  and  bottom  of  the  domain  we  specified  the  normal  derivative  of  u>.  This 
data  was  incorporated  into  the  difference  stendl  using  the  method  of  fictitious  points.  For 
points  on  such  boundaries,  the  central  difference  occurring  in  (2.9)  uses  a  point  just  outside 
the  computational  domain.  This  point  is  eliminated  by  using  the  normal  derivative  boundary 
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condition.  For  example,  at  points  on  the  plate,  the  second  derivative  term  is  approximated  by 

dy  ay  ^ 


dy 1 


For  the  data  at  the  bottom  we  used 


~-li.o  *  -7-  Si [  (2.12) 

i.e.  a  discrete  analog  of  (2.8).  At  the  top  we  used  the  boundary  condition  -  0-  This 

was  chosen  because  it  insures  that  the  boundary  condition  u*0aty*0is  satisfied  exactly 
by  the  discrete  scheme.  This  fact  can  be  verified  by  directly  calculating  the  difference  in  the 
velocity  at  the  plate  mHncwH  by  <iii+1  and  cd*.  The  procedure  is  a  discrete  version  of  the  argu¬ 
ments  used  to  find  die  conditions  that  the  rime  derivative  of  the  constraint  (2.7)  vanish.  An 
upper  boundary  condition  appears,  however,  because  the  computational  domain  is  finite.  The 
dependence  of  the  velocity  at  the  plate  upon  the  far  field  boundary  condition  illustrates  the 
point  that  one  must  be  careful  about  the  numerical  implementation  of  such  schemes  in 
unbounded  domains.  This  problem  with  far  field  boundary  conditions  is  also  discussed  within 
the  context  of  the  Navier  Stokes  equations  in  [1]. 

Our  computational  results  correspond  to  a  parameter  selection  of  v  =  .2,  rx  *  2.0, 
y„  =  4.0,  dx  »  dy  *  .1,  8f  =  .0125,  and  the  onset  velocity  U0  -  1.0.  The  values  of  param¬ 
eters  v,  U0,  and  y„  were  selected  to  insure  that  the  majority  of  the  vorticity  was  confined 
to  the  region  y  s  when  0s;si,,  The  timestep  was  chosen  so  that 

fir  <  min (‘^■»  No  instability  was  observed  for  this  choice  of  parameters. 

The  use  of  the  boundary  condition  (2.8)  that  the  initial  vorticity  induces  a  velo¬ 

city  field  which  satisfies  all  of  the  boundary  conditions.  The  initial  distribution  of  vorticity 
used, 
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o  ,  U° 
<  2* 


*  o 


i  —  0..m 
i  m  I.jh  0  <  j 


satisfies  this  assumption.  This  choice  of  initial  conditions  corresponds  to  impulsively  starting 
the  fluid  at  time  t=0  with  a  velocity  U0. 

The  development  of  the  boundary  layer  is  depicted  in  figures  2.1(a)  to  2.1(c).  At  time 
t  *  0  the  velocity  is  uniform  with  value  U0  except  at  the  surface  of  the  plate  where  it  is  zero. 
At  a  later  time,  r  -  2.5,  a  typical  boundary  layer  profile  is  emerging  near  the  leading  edge. 
Further  down  the  plate  the  velocity  field  is  still  uniformly  horizontal,  as  the  upstream  influ¬ 
ence  is  not  felt.  The  fluid  is  retarded,  however,  because  of  the  diffusion  of  vorticity  from  the 
wall.  At  a  much  later  time  t  =  10.0  the  solution  has  reached  steady  state  and  the  boundary 
layer  profile  is  depicted  in  figure  2.1(c). 

A  steady  state  solution  of  the  equations  is  the  Blasius  solution.  This  is  a  solution  which 

U  1 

is  self  similar  with  a  similarity  variable  t)  =  y(— )2,  i.e.  the  steady  state  solution  (u,v)  is 

xv 

given  by 


«(*o0  =  «4(t0  v(*o0  =  v4(ti) 


where  («4Cn),v4(n))  is  the  Blasius  solution.  A  comparison  of  the  steady  state  profile  and  the 
Blasius  solution  obtained  from  [10]  is  given  in  figure  2.2.  In  this  figure,  the  u  component  of 
the  computed  velocity  is  shown  using  arrows,  while  the  u  component  for  the  corresponding 
Blasius  solution  is  denoted  by  the  solid  line.  The  solution  is  plotted  for  positions  x=  .5,  x=1.0 
and  x*1.5  along  the  plate.  The  difference  between  the  computed  solution  and  the  Blasius 
solution  is  extremely  small. 

To  see  more  dearly  the  accuracy  with  which  the  Blasius  profile  is  approximated,  we 
show  in  figure  2.3  a  plot  of  the  difference  between  the  computed  «  velocity  and  the  u  velocity 
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of  the  Blasius  solution.  The  results  are  plotted  verses  the  similarity  variable  r\  so  that  each 
profile  can  be  compared  with  the  other.  The  velocity  profile  which  is  most  in  error  is  that 
near  the  leading  edge  of  the  plate,  x=.f .  This  is  to  be  expected  since  the  number  of  points 
per  profile  is  less  at  this  station  than  at  stations  further  down  the  plate.  However,  the  max¬ 
imum  relative  error  of  all  the  profiles  is  less  than  1%. 

These  results  demonstrate  that  it  is  possible  to  incorporate  the  boundary  condition  (2.8) 
into  a  numerical  scheme  and  obtain  accurate  answers.  Furthermore,  this  incorporation  can  be 
done  within  the  framework  of  an  explicit  method.  (There  has  been  some  doubt  about  the  pos¬ 
sibility  of  this.)  The  use  of  an  explicit  method  is  limiting  because  of  the  time  step  restriction. 
A  semi-implicit  finite  difference  scheme  has  been  developed  and  used  to  obtain  results  which 
agree  with  those  given  above.  We  chose  not  to  discuss  the  scheme  because  it  is  a  bit  more 
difficult  to  describe  and  not  any  more  illuminating  as  an  example  than  the  purely  explicit 
scheme. 

We  conclude  this  section  by  discussing  the  relation  between  a  method  introduced  by 
Chorin  far  solving  the  Prandtl  boundary  layer  equations  and  an  analogous  finite  difference 
scheme.  In  [4],  Chorin  uses  the  method  of  fractional  steps  and  vortidty  boundary  creation  to 
solve  equations  (2.1)-(2.5).  The  discretization  used  is  based  on  computational  elements  which 
are  segments  of  vortex  sheets.  The  basic  timestep,  ignoring  the  precise  implementation 
details,  are  as  follows:  An  approximate  solution  of  the  invisdd  equations  is  advanced  one 
time-step  (equation  2.1  without  v).  This  leads  to  a  vortidty  distribution  which  induces  a  velo 
dty  field  which  does  not  satisfy  the  tangential  boundary  condition  on  the  plate.  If  the  velodty 
field  at  location  idx  on  the  plate  is  U„  then  a  sheet  of  strength  2U,  is  created  on  the  boun¬ 
dary.  This  sheet  is  then  allowed  to  partidpate,  with  the  sheets  already  in  the  fluid,  in  a  ran¬ 
dom  walk.  This  random  walk  approximates  one  step  of  the  evolution  of  the  viscous  terms  of 
the  equation. 


Consider  the  following  fractional  step  scheme  which  uses  the  finite  difference  approxi¬ 


mations  given  previously, 
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(I) 
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LL  «  _ 


A,/**,*#4) 


"•J'-.ZU.  =  vD^D-^r1 


(2.13) 

(2.14) 


The  normal  derivative  boundary  condition 


du  .  - 1  £  r  A<.»  ♦  i(«**.  «*)  +  A^aSw4) , . 

171'-0  T,?„[ - 2 — - * 


(2.15) 


is  incorporated  into  the  finite  difference  stencil  occurring  in  (2.14).  The  upper  boundary  con¬ 


dition  is 


.  dd 


0.  One  can  show,  as  with  the  explicit  scheme  given  previously,  that  if  the 


dy 

velocity  field  induced  by  u*  satisfies  (2.5)  then  the  velocity  induced  by  w4+1  satisfies  (2.5)  as 
well. 


The  effect  of  the  boundary  condition  (2.15)  is  to  induce  vortirity  on  die  boundary  of  the 
computational  domain.  One  can  isolate  the  amount  of  vortidty  induced  by  this  boundary  con¬ 
dition  by  considering  one  step  of  equation  (2.14)  using  zero  initial  data  and  boundary  condi¬ 
tion  (2.15).  This  is  die  vortidty  u  defined  by 


<»ij  ~  «o, 


8  r 


L  ■  vd;d;<% 


with  initial  condition  =  0  and  boundary  condition  (2.15).  Working  through  the  algebra 
we  find  that  u  is  non-zero  only  at  the  lower  boundary,  and  is  given  by 


a) 


28r  £  r  Aiji  ♦  i(<*.  »*)  +  A^OJ*,®4) , . 

“  '  -  T& - 2 — - to 


(2.16) 


If  we  use  the  expression  for  A tJ  from  (2.13)  in  (2.16)  we  find  that 


,VtJ 


m 


S 


x(Uo +  2 

«y  ^  -o  2 


(2.17) 


Here  we  are  making  the  assumption  that  <nfj  satisfies  (2.5)  discretely,  i.e. 


Uo+  j. 

y-o  2 


The  quantity  U,  in  (2.17)  is  the  u  velocity  at  the  point  idx  on  the  plate  which  is  induced 


by  »*+1.  From  the  definition  of  w*+1  we  see  that  this  is  the  slip  induced  by  one  step  of  Euler 


2Ui 

flow.  The  vortidty  corresponds  to  a  computational  element  which  induces  a  jump  in 


vdodty  of  2U„  i.e.  a  sheet  with  strength  2U,.  This  is  precisely  the  strength  of  the  vortex 


sheet  which  is  introduced  by  Chorin  in  his  scheme.  Thus,  there  is  exact  agreement  in  the 


amount  of  vortidty  which  is  added  at  each  time  step.  We  believe  that  this  is  a  remarkable 


result  in  light  of  die  difference  in  the  types  of  discretizations  used. 


t  *  in*  *;» 


*  r 


BlbUography 


[1]  Anderson,  C.R.,  "Vortidty  Boundary  Conditions  and  Vortidty  Boundary  Generation  for 
Two  Dimensional  Incompressible  Flows” ,  in  preparation. 

[2]  Batchelor,  G,,  An  Introduction  to  Fluid  Mechanics.  Cambridge  University  Press,  Cam¬ 
bridge,  Great  Britain,  1967,  pg.  277. 

[3]  Chorin,  A.J.,  "Numerical  Study  of  Sightly  Viscous  Flow”,  J.  Fluid  Mech.,  57,  (1973), 
Pg-  785. 

[4]  Chorin.  A.J.,  "Vortex  Models  and  Boundary  Layer  Instability",  SLAM  J.  Sd.  Stat.  Corn- 
put.,  1.  1980),  pg.  1. 

[5]  ColeOa,  P,  "Multidimensional  Upwind  Schemes  for  Hyperbolic  Conservation  Laws",  sub¬ 
mitted  'or  publication,  J.  Comput.  Phys.,  1986. 

[6]  Gupta,  M.  and  Manohar,  R.,  "Boundary  Approximations  and  Accuracy  in  Viscous  Flow 
Computations",  J.  Comput.  Phys.,  31,  (1979),  pg.  265. 

[7]  Orszag,  S.  and  Israeli,  M.,  "Numerical  Simulation  of  Viscous  Incompressible  Flows”, 
Annual  Rev.  of  Fluid  Mech.,  6,  (1974),  pg.  281. 

[8]  QuartapeQe,  L.  and  Valz-Gris,  F.,  "Projection  Conditions  on  the  Vortidty  in  Viscous 
Incompressible  Flows",  Int.  J.  Numer.  Math.  Fluids,  1  (1981)  pg.  129. 

[9]  QuartapeDe,  L.,  "Vortidty  Conditioning  in  the  Computation  of  Two  Dimensional 
Viscous  Flows",  J.  Comput.  Phys.,  40.  (1981),  pg.  453. 

[10]  Schlkfating,  R,  Boundary  Layer  Theory.  McGraw-Hill,  New  York,  1960,  pg.  139. 


